
			/* This is program testxloop.cpp, to compare to onloop4pt.cpp, we write this program */
			/*       author: Phan Hong Khiem  
			/*		 Phan Le Anh Quan
			/*		 Do Hoang Son.                                                       */
			/*       email: phanhongkhiem@gmail.com                                              */
			/*										     */
			/******************************************************************************************
			**	Execution Plan									
			**			+-------------------------------+				
			**			|				|             
			**			|            Input P ,m
			**			|				|
			**			+-------------------------------+
			**					|
			**					V
			**			+-------------------------------+
			**			|				|
			**			|       calculate q{ij}	        |
			**			|                       	|
			**			|				|
			**			+-------------------------------+
			**                                      |
			**                                      V
			**			+-------------------------------+
			**			|	    step1		|
			**			| 				|
			**			|      a_{lk},b_{l}{k},c_{lk}   |
			**			|				|
			**			+-------------------------------+
			**                                      |
			**                                      V
			**			+-------------------------------+
			**			|	   step2		|
			**			|     A_{mlk}, B_{mlk},C_{mlk}  |
						|     ........F_{nmlk}......    |
			**			|				|
			**			+-------------------------------+
			**                                      |
			**                                      V
			**			+-------------------------------+
			**			|	   step 3		|
			**			| 				|
			**			|           D_{0}		|
			**			|				|
			**			+-------------------------------+
			**
***********************************************************************************************************************************************/
#include <iostream>
#include <getopt.h>
#include <sstream>
#include "ginac/ginac.h"
#include "ginac/flags.h"

using namespace std;
using namespace GiNaC;	

/* *****************************************************************************************************************************/
/*                            DECLARE FUNCTION FROM THETA TO LogAG FUNCTION                                                    */
/*******************************************************************************************************************************/
ex theta(ex );
ex sign(ex );
ex delta(ex);
ex Etafunction (ex, ex);
ex Rfunction(ex, ex);
ex LogACG(ex , ex , ex  ,ex );
ex LogARG(ex , ex , ex  ,ex );
ex LogAG(ex , ex , ex  ,ex );
ex ThetaG( ex , ex , ex , ex);
ex ThetaGC(ex , ex , ex , ex, ex);
void input_option(int argc, char *argv[], lst &p, lst &m);	

/*******************************************************************************************************************************/
					
int main(int argc, char *argv[])
{
/* *****************************************************************************************************************************/
/*                               DECLARE VARIABLE qij, ms, A[m][l][k] ... T4[n][m][l][k].                                      */
/*******************************************************************************************************************************/	

	ex rho=1e-30;
	int l,k,m,n;
	ex q10;
	ex q20,q21;
	ex q30,q31,q32;
//ex m1s,m2s,m3s,m4s;
	ex ms[5];
// declare a,b,c, d,AC,alpha  matrix.
	ex a[5][5];
	ex b[5][5];
	ex c[5][5];
	ex d[5][5];
	ex AC[5][5];
	ex alpha[5][5];
	ex f[5][5];
	ex fminus[5][5];
// declare A,B,C,D,beta,phi, Q,....f,g  matrix.
	ex A[5][5][5];
	ex B[5][5][5];
	ex C[5][5][5];
	ex D[5][5][5];
	ex beta[5][5][5];
	ex phi[5][5][5];

	ex Q[5][5][5];
	ex E[5][5][5];
	ex P[5][5][5];

	ex g[5][5][5];
	ex gminus[5][5][5];

	ex z1phi[5][5][5];
	ex z2phi[5][5][5];

	ex z1beta[5][5][5];
	ex z2beta[5][5][5];
	
	ex A0[5][5][5];
	ex B0[5][5][5];
	ex C0[5][5][5];
// declare F matrix, T1,2,3,4 ; Oplus, Ominus MATRIX
	ex F[5][5][5][5];
	ex T1[5][5][5][5];
	ex T2[5][5][5][5];

	ex Oplus[5][5][5][5];
	ex Ominus[5][5][5][5];
/* *****************************************************************************************************************************/
/*                                 CACULATE: INPUT PS AND MS[K]                                                                */
/*******************************************************************************************************************************/	

	//lst  p_,m_;
	//input_option(argc, argv, p_,m_);
	//ex Im;
	//for(Im=0.0;Im <=1000.0;Im+= 10.0){
	ex p1s =10.0;// p_.op(0);
	ex p2s =50.0;// p_.op(1);
	ex p3s = 10.0;//p_.op(2);
	ex p4s = 70.0;//p_.op(3);
	ex p12s =150.0;//p_.op(4);
	ex p23s = 10.0;//p_.op(5);
	
	ms[1] 	=  65610.0-0.0*I;//m_.op(0);
	ms[2]	=  82810.0-0.0*I;//-I*rho2;//m_.op(1);
	ms[3] 	=  65610.0-0.0*I;//-I*rho1;//m_.op(2);
        ms[4]	=  82810.0-0.0*I;//-I*rho2;//m_.op(3);
	
	ex m1s=ms[1],m2s=ms[2],m3s=ms[3],m4s=ms[4];
	//ms[4]=m4s;
/* *****************************************************************************************************************************/
/*                                CACULATE:  LORENT TRANFORMATION, INPUT TO XLOOP                                             */
/*******************************************************************************************************************************/
	ex  p12 = (p12s - p1s - p2s)/2;
	ex  p23 = (p23s - p2s - p3s)/2;
	ex  p13 = (p2s + p4s - p12s - p23s)/2;
	
	q10 = sqrt(p1s);
	q20 = p12/sqrt(p1s) + sqrt(p1s);
	q21 = sqrt( pow(p12/sqrt(p1s) + sqrt(p1s),2) - p12s );
	q30 = p13/sqrt(p1s) + p12/sqrt(p1s) + sqrt(p1s);
	q31 = (q30*q20 - p12s - p13 - p23)/q21;	
	q32 = sqrt( pow(q30,2) - pow(q31,2) - p4s );

	ex q11=0.0, q12=0.0, q13=0.0, q22=0.0, q23=0.0, q33=0.0,q40=0.0, q41=0.0, q42=0.0, q43=0.0;
/* *****************************************************************************************************************************/
/*                                 CACULATE:  a[l][K]    =2.0 (q[l][0]-q[k][0])                                                 */
/*******************************************************************************************************************************/
	a[1][1]=0.0;	
	a[1][2]=(q10-q20)*2.0;
	a[1][3]=(q10-q30)*2.0;
	a[1][4]=(q10-q40)*2.0;
	
/*for l=2, k change*/
	
	a[2][1]=(q20-q10)*2.0;
	a[2][3]=(q20-q30)*2.0;
	a[2][4]=(q20-q40)*2.0;
	a[2][2]=0.0;
/*for l=3, k change*/
	a[3][1]=(q30-q10)*2.0;
	a[3][2]=(q30-q20)*2.0;
	a[3][4]=(q30-q40)*2.0;
	a[3][3]=0.0;
/*for l=4, k change*/
	a[4][1]=(q40-q10)*2.0;
	a[4][2]=(q40-q20)*2.0;
	a[4][3]=(q40-q30)*2.0;
	a[4][4]=0.0;
/* *****************************************************************************************************************************/
/*                                 CACULATE:  b[l][K]    = -2.0*(q[l][1]-q[k][1])                                              */
/*******************************************************************************************************************************/
	
	b[1][1]=0.0;
	b[1][2]=-(q11-q21)*2.0;
	b[1][3]=-(q11-q31)*2.0;
	b[1][4]=-(q11-q41)*2.0;
	
	b[2][1]=-(q21-q11)*2.0;
	b[2][2]=0.0;
	b[2][3]=-(q21-q31)*2.0;
	b[2][4]=-(q21-q41)*2.0;
	
	b[3][1]=-(q31-q11)*2.0;
	b[3][2]=-(q31-q21)*2.0;
	b[3][3]=0.0;
	b[3][4]=-(q31-q41)*2.0;
	
	b[4][1]=-(q41-q11)*2.0;
	b[4][2]=-(q41-q21)*2.0;
	b[4][4]=0.0;
	b[4][3]=-(q41-q31)*2.0;
     
/* *****************************************************************************************************************************/
/*                                 CACULATE:  c[l][K]    = q[l][2]-q[k][2]                                                     */
/*******************************************************************************************************************************/
	c[1][1]=0;
	c[1][2]=-(q12-q22)*2.0;
	c[1][3]=-(q12-q32)*2.0;
	c[1][4]=-(q12-q42)*2.0;
	
/*for l=2, k c[hange*/
	
	c[2][1]=-(q22-q12)*2.0;
	c[2][2]=0.0;
	c[2][3]=-(q22-q32)*2.0;
	c[2][4]=-(q22-q42)*2.0;
	
/*for l=3, k c[hange */

	c[3][1]=-(q32-q12)*2.0;
	c[3][2]=-(q32-q22)*2.0;
	c[3][3]=0.0;
	c[3][4]=-(q32-q42)*2.0;

/*for l=4,k c[hange*/

	c[4][1]=-(q42-q12)*2.0;
	c[4][2]=-(q42-q22)*2.0;
	c[4][4]=0.0;
	c[4][3]=-(q42-q32)*2.0;

/* *****************************************************************************************************************************/
/*                                 CACULATE:  d[l][K]    = (q[l]-q[k])^2-m^2[k]-m^2[l]                                         */
/*******************************************************************************************************************************/
	
	d[1][1]=0;
	d[1][2]=pow(q10-q20,2)-pow(q21,2)-(m1s-m2s);
	d[1][3]=pow(q10-q30,2)-pow(q31,2)-pow(q32,2)-(m1s-m3s);
	d[1][4]=pow(q10,2)-(m1s-m4s);

	d[2][1]=pow(q20-q10,2)-pow(q21,2)-(m2s-m1s);
	d[2][2]=0.0;
	d[2][3]=pow(q20-q30,2) - pow(q21-q31,2) - pow(q32,2) - (m2s-m3s);
	d[2][4]=pow(q20,2)-pow(q21,2)-(m2s-m4s);
	
	d[3][1]=pow(q30-q10,2)-pow(q31,2)-pow(q32,2)-(m3s-m1s);
	d[3][2]=pow(q30-q20,2)-pow(q31-q21,2)-pow(q32,2)-(m3s-m2s);
	d[3][3]=0.0;
	d[3][4]=pow(q30,2)-pow(q31,2)-pow(q32,2)-(m3s-m4s);

	d[4][1]=pow(q10,2)-(m4s-m1s);
	d[4][2]=pow(q20,2)-pow(q21,2)-(m4s-m2s);
	d[4][4]=0.0;
	d[4][3]=pow(q30,2)-pow(q31,2)-pow(q32,2)-(m4s-m3s);
	
/* *****************************************************************************************************************************/
/*                                 CACULATE:  AC[l][K]    = a[l][k]+c[l][k]                                                    */
/*******************************************************************************************************************************/
 for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if(l!=k){
							AC[k][l]=a[k][l]+c[k][l];
						   } 
	
/* *****************************************************************************************************************************/
/*                                 CACULATE:  alpha[l][k]    = b[l][k]/AC[l][k]                                                */
/*******************************************************************************************************************************/
 for (k = 1; k<=4; k++) for(l = 1; l<=4; l++) if(l!=k){

					if(AC[k][l]==0){

						cout<<"alpha_"<< l << k << "have denorminator is zero (:????"<<endl;
					}	
					else{		
						alpha[k][l]=b[k][l]/AC[k][l];
					}
			} 
	
/* *****************************************************************************************************************************/
/*                                 CACULATE:  A[m][l][K]    = a[l][k]-a[l][k]*AC[m][k]/AC[l][k]                                */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
					if(AC[l][k]==0){
						cout<<"A_"<< m << l << k << "have denorminator is zero (:"<<endl;
					}	
					else{		
						A[m][l][k]=a[m][k]-a[l][k]*AC[m][k]/AC[l][k];
					}
			}   
	
/* *****************************************************************************************************************************/
/*                                 CACULATE: B[m][l][k]=b[m][k]-b[l][k]*AC[m][k]/AC[l][k];                                     */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(AC[l][k]==0){
					cout<<"B_"<< m << l << k << "have denorminator is zero (:"<<endl;
				}	
				else
				{		
					B[m][l][k]=b[m][k]-b[l][k]*AC[m][k]/AC[l][k];
				}
			}   
	
/* *****************************************************************************************************************************/
/*                                 CACULATE: C[m][l][k]=d[m][k]-d[l][k]*AC[m][k]/AC[l][k];                                      */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(AC[l][k]==0){
					cout<<"C_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else
				{		
					C[m][l][k]=d[m][k]-d[l][k]*AC[m][k]/AC[l][k];
				}
			}   
/* *****************************************************************************************************************************/
/*                                 CACULATE:     D[m][l][k]=1.0-2.0*a[l][k]/AC[l][k]+alpha[l][k]*alpha[l][k];  test....         */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(AC[l][k]==0){
					cout<<"D_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else
				{
					D[m][l][k]=1.0-2.0*a[l][k]/AC[l][k] + alpha[l][k]*alpha[l][k];
				}
	    		}

/* *****************************************************************************************************************************/
/*                                 CACULATE:   beta[m][l][k]=(
/*							(A[m][l][k]/B[m][l][k]-alpha[l][k])
/*						+sqrt((A[m][l][k]/B[m][l][k]-alpha[l][k])*(A[m][l][k]/B[m][l][k]-alpha[l][k])  */
/*							      -D[m][l][k] )						       */
/*						      )/D[m][l][k];                                                            */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(D[m][l][k]==0){
					cout<<"beta_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else {
					beta[m][l][k]=(
							(A[m][l][k]/B[m][l][k]-alpha[l][k])
							 +sqrt((A[m][l][k]/B[m][l][k]-alpha[l][k])*(A[m][l][k]/B[m][l][k]-alpha[l][k])
							       -D[m][l][k] )
						      )/D[m][l][k];
				}
}
/* *****************************************************************************************************************************/
/*                            CACULATE:     phi[m][l][k]=
					       (A[m][l][k]/B[m][l][k]-alpha[l][k])
					       +sqrt((A[m][l][k]/B[m][l][k]-alpha[l][k])*(A[m][l][k]/B[m][l][k]-alpha[l][k])
					      -D[m][l][k] );                                                                   */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){

				phi[m][l][k]=( A[m][l][k]/B[m][l][k]-alpha[l][k])
					        +sqrt((A[m][l][k]/B[m][l][k]-alpha[l][k])*(A[m][l][k]/B[m][l][k]-alpha[l][k])
					        -D[m][l][k] );
}
	
/* *****************************************************************************************************************************/
/*                            CACULATE:  Q[m][l][k]=-2.0*(C[m][l][k]/B[m][l][k]+d[l][k]*beta[m][l][k]/AC[l][k]);               */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(AC[l][k]==0 || B[m][l][k]==0){
					cout<<"Q_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else {
					Q[m][l][k]=-2.0*(C[m][l][k]/B[m][l][k]+d[l][k]*beta[m][l][k]/AC[l][k]);
				}
			}
/* *****************************************************************************************************************************/
/*                            CACULATE:  E[m][l][k]=-2.0*(C[m][l][k]*phi[m][l][k]/B[m][l][k]+d[l][k]/AC[l][k] );               */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if(AC[l][k]==0 || B[m][l][k]==0){
					cout<<"E_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else {
					E[m][l][k]=-2.0*(C[m][l][k]*phi[m][l][k]/B[m][l][k]+d[l][k]/AC[l][k] );
				}
			}
/* *****************************************************************************************************************************/
/*                            CACULATE: P[m][l][k]=-2.0*( 
							  (A[m][l][k]/B[m][l][k]-alpha[l][k])*(1+beta[m][l][k]*phi[m][l][k])
							  -D[m][l][k]*beta[m][l][k]-phi[m][l][k]
							);            							       */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k){
				if( B[m][l][k]==0){
					cout<<"P_"<< m << l << k << "have denorminator is zero (:????"<<endl;
				}	
				else{
					P[m][l][k]=-2.0*( 
							  (A[m][l][k]/B[m][l][k]-alpha[l][k])*(1+beta[m][l][k]*phi[m][l][k])
							  -D[m][l][k]*beta[m][l][k]-phi[m][l][k]
							);
				}
}
	
/* *****************************************************************************************************************************/
/*                         CACULATE:  f[lk]				                                                      */
/*******************************************************************************************************************************/

for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k){
		
				ex temp=imag_part(-d[l][k]/AC[l][k]);	
		
				if (temp>0.0) {
					f[l][k] = 2.0;
				}			
				if (temp==0.0) {
					f[l][k] = 1.0;
				}	
				if (temp <0.0 ){
					f[l][k] = 0.0;
				}
			}

/* *****************************************************************************************************************************/
/*                         CACULATE:  fminus[l][k]			                                                      */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k)for(m = 1; m<=4; m++) if(m!=l && m!=k){
		
				ex temp=imag_part(-d[l][k]/AC[l][k]);	
		
				if (temp<0) {
					fminus[l][k] = 2.0;
				}			
				if (temp==0) {
					fminus[l][k] = 1.0;
				}	
				if (temp >0 ){
					fminus[l][k] = 0.0;
				}
			}

/* *****************************************************************************************************************************/
/*                         CACULATE:  g[mlk]                                                 					 */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k)for(m = 1; m<=4; m++) if(m!=l && m!=k){
		
				ex temp1=imag_part(-C[m][l][k]/B[m][l][k]);	
		
				if (temp1 > 0) {
					g[m][l][k] = 2.0;
				}			
				if (temp1==0) {
					g[m][l][k] =  1.0;
				}	
				if (temp1 <0 ){
					g[m][l][k] = 0.0;
				}
		
			}

/* *****************************************************************************************************************************/
/*                         CACULATE:  gminus[mlk]                                                                                  */
/*******************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k)for(m = 1; m<=4; m++) if(m!=l && m!=k){
		
				ex temp1=imag_part(-C[m][l][k]/B[m][l][k]);	
		
				if (temp1<0) {
					gminus[m][l][k] = 2.0;
				}			
				if (temp1==0) {
					gminus[m][l][k] = 1.0;
				}	
				if (temp1 >0 ){
					gminus[m][l][k] = 0.0;
				}
		
			}
/* *******************************************************************************************************************************************/
/*                            CACULATE: z1beta[m][l][k]
				   	                                                                                                     */
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {

		if(P[m][l][k]*beta[m][l][k]==0){

			cout<<"z1beta_"<< m << l << k << "have denorminator is zero (:"<<endl;

		}	
		else{
		     z1beta[m][l][k]=(
		        -(E[m][l][k]-Q[m][l][k]/beta[m][l][k])
			+sqrt( (E[m][l][k]-Q[m][l][k]/beta[m][l][k])*(E[m][l][k]-Q[m][l][k]/beta[m][l][k])-4.0*P[m][l][k]*(ms[k]-I*rho)/beta[m][l][k])
					)/(-2.0*P[m][l][k]/beta[m][l][k]);
			}
}

/* *******************************************************************************************************************************************/
/*                            CACULATE: z2beta[m][l][k]                                                                                      */
/*********************************************************************************************************************************************/
 for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {

				if(P[m][l][k]*beta[m][l][k]==0){

					cout<<"z2beta_"<< m << l << k << "have denorminator is zero (:"<<endl;

				}	
				else{
				     z2beta[m][l][k]=(
				        -(E[m][l][k]-Q[m][l][k]/beta[m][l][k])
					-sqrt( (E[m][l][k]-Q[m][l][k]/beta[m][l][k])*(E[m][l][k]-Q[m][l][k]/beta[m][l][k])
					      -4.0*P[m][l][k]*(ms[k]-I*rho)/beta[m][l][k])
							)/(-2.0*P[m][l][k]/beta[m][l][k]);
					}
}

/***********************************************************************************************************************/
 /*                             calculate z1phi [mlk] and z2phi[mlk]                                               *   /***********************************************************************************************************************/											                     

for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {
		if(P[m][l][k]*phi[m][l][k]==0){
			cout<<"z1phi_"<< m << l << k << "have denorminator is zero (:"<<endl;
		}	
		else{
		     z1phi[m][l][k]=(
			-(E[m][l][k]-Q[m][l][k]*phi[m][l][k])
			+sqrt( (E[m][l][k]-Q[m][l][k]*phi[m][l][k])*(E[m][l][k]-Q[m][l][k]*phi[m][l][k])-4.0*P[m][l][k]*(ms[k]-I*rho)*phi[m][l][k])
				    )/(-2.0*P[m][l][k]*phi[m][l][k]);
		}
}
/* *******************************************************************************************************************************************/
/*                            CACULATE:    z2phi[m][l][k]
               														                     */
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {

		if(P[m][l][k]*phi[m][l][k]==0){

			cout<<"z2phi_"<< m << l << k << "have denorminator is zero (:"<<endl;

		}	
		else{

		     z2phi[m][l][k]=(
			 -(E[m][l][k]-Q[m][l][k]*phi[m][l][k])
			 -sqrt( (E[m][l][k]-Q[m][l][k]*phi[m][l][k])*(E[m][l][k]-Q[m][l][k]*phi[m][l][k])-4.0*P[m][l][k]*(ms[k]-I*rho)*phi[m][l][k])
				    )/(-2.0*P[m][l][k]*phi[m][l][k]);

		}
}
/* *******************************************************************************************************************************************/
/*                            CACULATE:    A0[m][l][k]
               														                     */
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {
    A0[m][l][k]= imag_part(P[m][l][k]*E[m][l][k]);
}
/* *******************************************************************************************************************************************/
/*                            CACULATE:    B0[m][l][k]
               														                     */
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {
    B0[m][l][k]= imag_part(Q[m][l][k].conjugate()*E[m][l][k]-ms[k]*P[m][l][k]+I*rho*P[m][l][k]);
}
/* *******************************************************************************************************************************************/
/*                            CACULATE:   C0[m][l][k]
               														                     */
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k) {
    C0[m][l][k]= imag_part(-Q[m][l][k].conjugate()*ms[k]+Q[m][l][k].conjugate()*I*rho);
}

/* *******************************************************************************************************************************************/
/*           CACULATE: F[n][m][l][k]=(C[n][l][k]*B[m][l][k]-C[m][l][k]*B[n][l][k])/(A[n][l][k]*B[m][l][k]-A[m][l][k]*B[n][l][k]);            */ 											                      
/*********************************************************************************************************************************************/
	
		for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
			if(A[n][l][k]*B[m][l][k]-A[m][l][k]*B[n][l][k]==0){
				cout<<"F_"<<n<< m << l << k << "have denorminator is zero (:????"<<endl;
			}	
			else
			{
			 F[n][m][l][k]=(C[n][l][k]*B[m][l][k]-C[m][l][k]*B[n][l][k])/(A[n][l][k]*B[m][l][k]-A[m][l][k]*B[n][l][k])+I*rho*beta[m][l][k];
			}
		   }
/* *******************************************************************************************************************************************/
/*                                CACULATE: THE VALUE OF T1[n][m][l][k]                                                                      */ 																					                      
/*********************************************************************************************************************************************/
	for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
		if(P[m][l][k]==0.0){
			cout<<"T1_"<<n<< m << l << k << "have denorminator is zero (:????"<<endl;
		}	
		else
		{
		 T1[n][m][l][k]=( (Q[m][l][k]+P[m][l][k]*F[n][m][l][k]-beta[m][l][k]*E[m][l][k])+
			        sqrt( pow((Q[m][l][k]+P[m][l][k]*F[n][m][l][k]-beta[m][l][k]*E[m][l][k]),2.0)
				-4.0*P[m][l][k]*(Q[m][l][k]*F[n][m][l][k] +beta[m][l][k]*ms[k]-I*beta[m][l][k]*rho  ) ) )/(-2.0*P[m][l][k]);
		}
	   }
/* *******************************************************************************************************************************************/
/*                                CACULATE: THE VALUE OF T2[n][m][l][k]                                                                      */ 																					                      
/*********************************************************************************************************************************************/

	for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
		if(P[m][l][k]==0.0){
			cout<<"T2_"<< n << m << l << k << "have denorminator is zero (:????"<<endl;
		}	
		else
		{
		T2[n][m][l][k]=( (Q[m][l][k]+P[m][l][k]*F[n][m][l][k]-beta[m][l][k]*E[m][l][k])-
			       sqrt( pow((Q[m][l][k]+P[m][l][k]*F[n][m][l][k]-beta[m][l][k]*E[m][l][k]),2.0)
			       -4.0*P[m][l][k]*(Q[m][l][k]*F[n][m][l][k] +beta[m][l][k]*ms[k]-I*beta[m][l][k]*rho ) ))/(-2.0*P[m][l][k]);
		}
	   }


/* *******************************************************************************************************************************************/
/*                                CACULATE: THE VALUE OF Oplus[n][m][l][k]                                                                   */ 																					                      
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
		if(beta[m][l][k]==0.0){
			cout<<"Oplus_"<<n<< m << l << k << "have denorminator is zero (:????"<<endl;
		}	
		else
		{
		 Oplus[n][m][l][k]=( f[l][k]*g[m][l][k] + fminus[l][k]*g[m][l][k] )*log((F[n][m][l][k])/beta[m][l][k])

		   -( f[l][k]*g[m][l][k] + fminus[l][k]*g[m][l][k] )*theta( imag_part(-P[m][l][k]*z1beta[m][l][k]/beta[m][l][k]) )
						    *theta( imag_part( z2beta[m][l][k] ) )*2.0*Pi*I
		   +f[l][k]*g[m][l][k]*theta( imag_part(-P[m][l][k]*z1phi[m][l][k]*phi[m][l][k]) )
				      *theta( imag_part( z2phi[m][l][k] ) )*2.0*Pi*I
		   -f[l][k]*gminus[m][l][k]*theta( imag_part(-P[m][l][k]*z1phi[m][l][k]*phi[m][l][k]) )
				      *theta( - imag_part( z2phi[m][l][k] ) )*2.0*Pi*I
	    	 	;
		}
}
/* *******************************************************************************************************************************************/
/*                                CACULATE: THE VALUE OF Ominus[n][m][l][k]                                                                      */ 																					                      
/*********************************************************************************************************************************************/
for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
		if(beta[m][l][k]==0.0){
			cout<<"Ominus_"<<n<< m << l << k << "have denorminator is zero (:????"<<endl;
		}	
		else
		{
		Ominus[n][m][l][k] =  -fminus[l][k]*gminus[m][l][k]*log((F[n][m][l][k])/beta[m][l][k])
				      -f[l][k]*gminus[m][l][k]*log(-(F[n][m][l][k])/beta[m][l][k])
				      +fminus[l][k]*gminus[m][l][k]*theta( imag_part(-P[m][l][k]*z1beta[m][l][k]/beta[m][l][k]) )
						    *theta( imag_part( z2beta[m][l][k] ) )*2.0*Pi*I
				      -f[l][k]*gminus[m][l][k]*theta( imag_part(-P[m][l][k]*z1beta[m][l][k]/beta[m][l][k]) )
						    *theta(- imag_part( z2beta[m][l][k] ) )*2.0*Pi*I
				      -(fminus[l][k]*gminus[m][l][k]+fminus[l][k]*g[m][l][k])*theta( imag_part(-P[m][l][k]*z1phi[m][l][k]*phi[m][l][k]) )
				      							  *theta( imag_part( z2phi[m][l][k] ) )*2.0*Pi*I	;
		}
 }

// THE MAIN CALCULATE HERE.

ex D0=0.0;
ex term1=0.0;
ex term2=0.0;
ex term3=0.0;
ex term4=0.0;

for(k = 1; k<=4; k++) for(l = 1; l<=4; l++) if (l!=k) for(m = 1; m<=4; m++) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){

ex jacobian = abs(1.0-beta[m][l][k]*phi[m][l][k])*(1.0-delta(AC[l][k]))*(1.0-delta(B[m][l][k]))*(1.0/AC[l][k])
	          *(1.0/(B[m][l][k]*A[n][l][k]-B[n][l][k]*A[m][l][k]))*(-1.0/P[m][l][k]);

//term1 la so hang khong chua log.
term1= term1 + (Oplus[n][m][l][k]*Rfunction(-T1[n][m][l][k],-T2[n][m][l][k])+Ominus[n][m][l][k]*Rfunction(T1[n][m][l][k],T2[n][m][l][k]) )*jacobian;

// term2 la2 so hang chua F[nmlk]
term2= term2 
	    +(-f[l][k]*g[m][l][k]*LogAG( (1.0-beta[m][l][k]*phi[m][l][k])/beta[m][l][k],( F[n][m][l][k])/beta[m][l][k], -T1[n][m][l][k],-T2[n][m][l][k] )
            +(fminus[l][k]*gminus[m][l][k]+fminus[l][k]*g[m][l][k])*LogAG(-(1.0-beta[m][l][k]*phi[m][l][k])/beta[m][l][k],(F[n][m][l][k])/beta[m][l][k],T1[n][m][l][k],T2[n][m][l][k])  )*jacobian  
 	    -f[l][k]*gminus[m][l][k]*LogAG( -(1.0-beta[m][l][k]*phi[m][l][k])/beta[m][l][k], -(F[n][m][l][k])/beta[m][l][k], -T1[n][m][l][k],-T2[n][m][l][k] )*jacobian;

// term 3 la so hang chua log nhung ko co F[nmlk]
 term3= term3+
	  (     -(f[l][k]*g[m][l][k]+fminus[l][k]*g[m][l][k])*LogAG(-P[m][l][k]/beta[m][l][k],P[m][l][k]*z1beta[m][l][k]/beta[m][l][k], -T1[n][m][l][k],-T2[n][m][l][k] )


           -(f[l][k]*g[m][l][k]+fminus[l][k]*g[m][l][k])*LogAG(1.0,-z2beta[m][l][k], -T1[n][m][l][k],-T2[n][m][l][k] )


           + f[l][k]*g[m][l][k]*LogAG(-P[m][l][k]*phi[m][l][k], P[m][l][k]*phi[m][l][k]*z1phi[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] ) 


	   + f[l][k]*g[m][l][k]*LogAG(1.0, -z2phi[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] ) 


	   + f[l][k]*gminus[m][l][k]*LogAG(P[m][l][k]*phi[m][l][k], -P[m][l][k]*phi[m][l][k]*z1phi[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] )


           + f[l][k]*gminus[m][l][k]*LogAG(1.0, -z2phi[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] )


	   +(fminus[l][k]*g[m][l][k]-f[l][k]*gminus[m][l][k])*LogAG(P[m][l][k],Q[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] ) 
		//so hang co can -infty to 0

            + fminus[l][k]*gminus[m][l][k]*LogAG( P[m][l][k]/beta[m][l][k], P[m][l][k]*z1beta[m][l][k]/beta[m][l][k], T1[n][m][l][k],T2[n][m][l][k] ) 


	    +fminus[l][k]*gminus[m][l][k]*LogAG( -1.0, -z2beta[m][l][k], T1[n][m][l][k],T2[n][m][l][k] ) 


	    +f[l][k]*gminus[m][l][k]*LogAG( -P[m][l][k]/beta[m][l][k], -P[m][l][k]*z1beta[m][l][k]/beta[m][l][k], T1[n][m][l][k],T2[n][m][l][k] ) 


	    +f[l][k]*gminus[m][l][k]*LogAG( -1.0, -z2beta[m][l][k], T1[n][m][l][k],T2[n][m][l][k] ) 


	    -(fminus[l][k]*gminus[m][l][k]+fminus[l][k]*g[m][l][k] )* LogAG( P[m][l][k]*phi[m][l][k], P[m][l][k]*phi[m][l][k]*z1phi[m][l][k],T1[n][m][l][k],T2[n][m][l][k] ) 


	    -(fminus[l][k]*gminus[m][l][k]+fminus[l][k]*g[m][l][k] )* LogAG(-1.0, -z2phi[m][l][k],T1[n][m][l][k],T2[n][m][l][k] ) 

	     +(fminus[l][k]*g[m][l][k]-f[l][k]*gminus[m][l][k])*LogAG(-P[m][l][k],Q[m][l][k],T1[n][m][l][k],T2[n][m][l][k] )         )*jacobian;
           // for m real (so hang extra term ) dong gop cua tich phan \eta[S/pz+Q]G(z)
//term4 = term4 +ThetaG( -P[m][l][k], -Q[m][l][k], -T1[n][m][l][k],-T2[n][m][l][k] )*2.0*Pi*I*jacobian;

           // for m complex.

term4= term4  + ( theta(-imag_part(Q[m][l][k]))* fminus[l][k]*g[m][l][k]  + theta(imag_part(Q[m][l][k]))* f[l][k]*gminus[m][l][k] )*
               ThetaGC(-A0[m][l][k],-B0[m][l][k], -C0[m][l][k],-T1[n][m][l][k],-T2[n][m][l][k] )*2.0*Pi*I*jacobian;
        }
D0=term1 + term2 + term3 + term4;
/*cout<< term1.evalf()<<endl;
cout<<term2.evalf()<<endl;
cout<<term3.evalf()<<endl;
cout<<term4.evalf()<<endl;*/
// COUT TO TEST

cout<<real_part(D0).evalf()<<"\t"<<imag_part(D0).evalf()<<endl;
//cout <<"("<<real_part(D0).evalf()<<"," <<imag_part(D0).evalf()<<")"<<endl;
//}
cout << "Khiem's Code" << endl << "q_ij" << endl;
		cout << q10 << "\t" << q11 << "\t" << q12 << "\t" << q13 << "\t" << endl;
		cout << q20 << "\t" << q21 << "\t" << q22 << "\t" << q23 << "\t" << endl;
		cout << q30 << "\t" << q31 << "\t" << q32 << "\t" << q33 << "\t" << endl;
		cout << q40 << "\t" << q41 << "\t" << q42 << "\t" << q43 << "\t" << endl;
		cout << "m_i square" << endl;
		cout << m1s << "\t" << m2s << "\t" << m3s << "\t" << m4s << "\t" << endl;
	
		cout << "terms\t real part \t image part" << endl;
		for(k=1; k<=4; k++){
			for(l=1; l<=4; l++) if(l!=k) {
				cout << "a_" << k << l << "\t" << real_part(a[k][l]) << "\t" << imag_part(a[k][l]) << endl;
			}
		}	
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "b_" << k << l << "\t" << real_part(b[k][l]) << "\t" << imag_part(b[k][l]) << endl;
		}		 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "c_" << k << l << "\t" << real_part(c[k][l]) << "\t" << imag_part(c[k][l]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "d_" << k << l << "\t" << real_part(d[k][l]) << "\t" << imag_part(d[k][l]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "AC_" << k << l << "\t" << real_part(AC[k][l]) << "\t" << imag_part(AC[k][l]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "alpha_" << k << l << "\t" << real_part(alpha[k][l]) << "\t" << imag_part(alpha[k][l]) << endl;
		}
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "f_" << k << l << "\t" << real_part(f[k][l]) << "\t" << imag_part(f[k][l]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) {
			cout << "fminus_" << k << l << "\t" << real_part(fminus[k][l]) << "\t" << imag_part(fminus[k][l]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "A_" << k << l << m << "\t" << real_part(A[k][l][m]) << "\t" << imag_part(A[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "B_" << k << l << m << "\t" << real_part(B[k][l][m]) << "\t" << imag_part(B[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "C_" << k << l << m << "\t" << real_part(C[k][l][m]) << "\t" << imag_part(C[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "D_" << k << l << m << "\t" << real_part(D[k][l][m]) << "\t" << imag_part(D[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "beta_" << k << l << m << "\t" << real_part(beta[k][l][m]) << "\t" << imag_part(beta[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "phi_" << k << l << m << "\t" << real_part(phi[k][l][m]) << "\t" << imag_part(phi[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "Q_" << k << l << m << "\t" << real_part(Q[k][l][m]) << "\t" << imag_part(Q[k][l][m]) << endl;
		}
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "P_" << k << l << m << "\t" << real_part(P[k][l][m]) << "\t" << imag_part(P[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "E_" << k << l << m << "\t" << real_part(E[k][l][m]) << "\t" << imag_part(E[k][l][m]) << endl;
		}
		for(k=1; k<=4; k++) for(l=1; l<=4; l++)if(l!=k) for(m=1; m<=4; m++)if(m!=l && m!=k) {
			cout << "g_" << k << l << m << "\t" << real_part(g[k][l][m]) << "\t" << imag_part(g[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "gminus_" << k << l << m << "\t" << real_part(gminus[k][l][m]) << "\t" << imag_part(gminus[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if (l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) for(n=1; n<=4; n++) if(n!=l && n!=k && n!=m){
			cout << "F_" << k << l << m << n << "\t" << real_part(F[k][l][m][n]) << "\t" << imag_part(F[k][l][m][n]) << endl;
		}
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "z1beta_" << k << l << m << "\t" << real_part(z1beta[k][l][m]) << "\t" << imag_part(z1beta[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "z2beta_" << k << l << m << "\t" << real_part(z2beta[k][l][m]) << "\t" << imag_part(z2beta[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "z1phi_" << k << l << m << "\t" << real_part(z1phi[k][l][m]) << "\t" << imag_part(z1phi[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if(l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) {
			cout << "z2phi_" << k << l << m << "\t" << real_part(z2phi[k][l][m]) << "\t" << imag_part(z2phi[k][l][m]) << endl;
		} 
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if (l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) for(n=1; n<=4; n++) if(n!=l && n!=k && n!=m){
			cout << "T1_" << k << l << m << n << "\t" << real_part(T1[k][l][m][n]) << "\t" << imag_part(T1[k][l][m][n]) << endl;
		}
		for(k=1; k<=4; k++) for(l=1; l<=4; l++) if (l!=k) for(m=1; m<=4; m++) if(m!=l && m!=k) for(n=1; n<=4; n++) if(n!=l && n!=k && n!=m){
			cout << "T2_" << k << l << m << n << "\t" << real_part(T2[k][l][m][n]) << "\t" << imag_part(T2[k][l][m][n]) << endl;
		}
		for(k=1; k<=4; ++k) for(l=1; l<=4; ++l) if (l!=k) for(m=1; m<=4; ++m) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
	        cout << "Ominus["<< k << l << m << n <<"]"<< "\t" << real_part(Ominus[k][l][m][n])<< "\t" << imag_part(Ominus[k][l][m][n]).evalf() << endl;
		}
		for(k=1; k<=4; ++k) for(l=1; l<=4; ++l) if (l!=k) for(m=1; m<=4; ++m) if(m!=l && m!=k)for(n= 1; n<=4; n++) if(n!=l && n!=k && n!=m){
	        cout << "Oplus["<< k << l << m << n <<"]"<< "\t" << real_part(Oplus[k][l][m][n])<< "\t" << imag_part(Oplus[k][l][m][n]).evalf() << endl;
		}
		
return 0;

}
/* *****************************************************************************************************************************/
/*                            theta function defind as theta(x)= 1 if x>=0; and =0 if x<0;                                     */
/*******************************************************************************************************************************/
ex theta(ex x){

	ex temp;

	if(x>=0.0){
        	temp=1.0;
	}
	else {
		temp=0.0;
	}
return temp;

}
ex sign(ex x){

	ex temp;

	if(x>=0.0){
        	temp=1.0;
	}
	else {
		temp=-1.0;
	}
return temp;

}


/* *****************************************************************************************************************************/
/*                            CACULATE:   delta function defind as delta(x)=0 if x<>0; and =1 if x==0.                         */
/*******************************************************************************************************************************/

ex delta(ex x){
	ex temp;
	
	if (x!=0.0){
		    temp=0.0;
	}
	else {
		    temp=1.0;
	}

 return temp;
}
						
/* *****************************************************************************************************************************/
/*                            CACULATE:   define Etafunction here                					       */
/*******************************************************************************************************************************/
ex Etafunction (ex x, ex y){
    ex temp;
    if(imag_part(x)==0.0&& imag_part(y)!=0.0){
                                             if(real_part(x)<0.0){
                                               temp =-2.0*Pi*I*theta(imag_part(y));}
					     else{
                                               temp=0.0;
                                                          }
				}
    if(imag_part(y)==0.0&& imag_part(x)!=0.0){
                                             if(real_part(y)<0.0){
                                               temp =-2.0*Pi*I*theta(imag_part(x));}
					     else{
                                               temp=0.0;
                                                          }
				}
     if(imag_part(y)!=0.0&& imag_part(x)!=0.0){
    temp=(  theta(-imag_part(x)) *theta(-imag_part(y))*theta(imag_part(x*y))
	     - theta(imag_part(x)) *theta(imag_part(y))*theta(-imag_part(x*y)) )*2.0*I*Pi;
                     }
   return temp;
}		
/* *****************************************************************************************************************************/
/*                            CACULATE:    define Rfunction here.           					       		*/
/*******************************************************************************************************************************/
ex Rfunction(ex x, ex y){
	ex temp;
	if(x==y){
		 cout<<"R function have zero's denominate "<< endl;
	}
	else{
	         temp =  (log(x)-log(y))/(x-y);
	}

   return temp;
}

/* *****************************************************************************************************************************/
/*                            CACULATE:    define    LogACG(ex a, ex b, ex x1 ,ex y)    			               */
/*******************************************************************************************************************************/
/*ex LogACG(ex a, ex b, ex x1 ,ex y){
	ex t=b/a;
	ex A=sqrt(real_part(x1)*real_part(x1)+imag_part(x1)*imag_part(x1)  )
	     +sqrt(real_part(y)*real_part(y)+imag_part(y)*imag_part(y)  )
	     +sqrt(real_part(t)*real_part(t)+imag_part(t)*imag_part(t)  );

	ex x0=x1/A;
	ex y0=y/A;
	ex z0=t/A;

	ex GiNaCLogACG=(  (log(a*A) )*(log(x0)-log(y0))/(x0-y0)
			+( -log(x0)*log(x0)/2.0 +log(y0)*log(y0)/2.0 +Li2(1.0-z0/y0)-Li2(1.0-z0/x0) )/(-x0+y0)
			+log(y0)*( Etafunction(z0-y0, 1.0/(1.0-y0)) - Etafunction(z0-y0, 1.0/(-y0))   )/(-x0+y0)
			-log(x0)*( Etafunction(z0-x0, 1.0/(1.0-y0)) - Etafunction(z0-x0, 1.0/(-x0))   )/(-x0+y0)
			+log(1.0-z0/y0)*Etafunction(z0,1.0/y0)/(-x0+y0)
			-log(1.0-z0/x0)*Etafunction(z0,1.0/x0)/(-x0+y0)  )/(A);

	return GiNaCLogACG;

}*/
/* *****************************************************************************************************************************/
/*                            CACULATE:    define    LogARG(ex a, ex b, ex x1 ,ex y)    			               */
/*******************************************************************************************************************************/
/*ex LogARG(ex a, ex b, ex x1 ,ex y){

	ex t=b/a;
	ex A=sqrt(real_part(x1)*real_part(x1)+imag_part(x1)*imag_part(x1)  )
	     +sqrt(real_part(y)*real_part(y)+imag_part(y)*imag_part(y)  )+
	      sqrt(real_part(t)*real_part(t)+imag_part(t)*imag_part(t)  );

	ex x0=x1/A;
	ex y0=y/A;
	ex z0=t/A;
	ex GiNaCLogARG=log(b)*(log(x1)-log(y))/(x1-y)+LogACG(a/b, 1.0,  x1, y);

	return GiNaCLogARG;

}*/
/* *****************************************************************************************************************************/
/*                            CACULATE:    define    LogAG(ex a, ex b, ex x1 ,ex y)    			                       */
/*******************************************************************************************************************************/
/*ex LogAG(ex a, ex b, ex x1 ,ex y){
			
	ex temp;
	
	if (a>0.0){		
		   temp =	LogACG(a, b, x1 , y);
	}
	else {
		   temp =	LogARG( a, b,  x1 , y);
	}

	return temp;
}*/
/*ex LogAG(ex a, ex b, ex x0 ,ex y0){
			
	ex temp;
	 ex z0=b/a;
	 temp=log(  a+I*1e-10*sign( imag_part(b) ))*Rfunction( x0, y0)+( log(x0)*log( (z0-x0)/(1.0-x0)  )-log(z0)*log ( (x0-z0)/x0 )+Li2(1.0/x0)
		-Li2(1.0-x0)-Li2(z0/x0) -log(y0)*log( (z0-y0)/(1.0-y0) )+log(z0)*log( (y0-z0)/(y0) )-Li2(1.0/y0)+Li2(1.0-y0) +Li2(z0/y0)  )/(x0-y0); 

	return temp;
}*/
ex LogAG(ex a , ex b, ex T1 ,ex T2){
		ex r=a/b;
		ex temp=( Li2(1.0-r*T1)-Li2(1.0-r*T2)+Etafunction(T1,r)*log(1.0-r*T1) - Etafunction(T2,r)*log(1.0-r*T2) )/(-T1+T2)
			+(log(T1)-log(T2))*log(b)/(T1-T2);
		return temp;
	}



ex ThetaG(ex a, ex b, ex x ,ex y){
ex temp;
  if(a==0){
		if(b>=0){
                	temp= Rfunction(x,y)+ Rfunction(-x,-y);
		}
	
               else{
		temp=0.0;
		}
	}

  if (a>0){
		temp= Rfunction(-b/a+x,-b/a+y);
		}
  if (a<0){ 
		temp = Rfunction(b/a-x,b/a-y);
	}
   return temp;

}
ex ThetaGC(ex a, ex b, ex c, ex x, ex y){

ex temp;
ex delta= b*b-4.0*a*c;
if(a==0.0&&b==0.0){
                 if(c>=0.0){
			    temp= Rfunction(x,y)+ Rfunction(-x,-y);
			      }
		 else  { 
			temp=0.0;}
  }
 if(b>0.0){   
              temp= Rfunction(-c/b+x,-c/b+y);
                 }        
 if(b<0.0){   
              temp= Rfunction(c/b-x,c/b-y);
                 }      
 if(a>0.0&& delta<=0.0){
	       temp= Rfunction(x,y)+ Rfunction(-x,-y);
                               }
if(a>0.0&& delta>0.0){
	ex z01=(-b+sqrt(delta) )/(2.0*a);
	ex z02=(-b-sqrt(delta) )/(2.0*a);
        temp= Rfunction(z01+x,z01+y)+ Rfunction(-z02-x,-z02-y);
}

if(a<0.0&& delta<=0){
	       temp= 0.0;
                               }
if(a<0.0&& delta>0.0){
	ex z01=(-b+sqrt(delta) )/(2.0*a);
	ex z02=(-b-sqrt(delta) )/(2.0*a);
        temp= Rfunction(z02+x,z02+y)+ Rfunction(z01+x,z01+y);
}
          
return temp;
}

/* *****************************************************************************************************************************/
/*                            CACULATE:    input_option(int argc, char *argv[], lst &p, lst &m)			               */
/*******************************************************************************************************************************/
void input_option(int argc, char *argv[], lst &p, lst &m){
		//Use optget_long to get the long options with their values
	ex m1, m2, m3, m4;
	m1 = symbol("m1"); m2 = symbol("m2"); m3 = symbol("m3"); m4 = symbol("m4");
	
	// Store the P's for loop tool
	ex p1s=symbol("p1s"),p2s=symbol("p2s"),p3s=symbol("p3s"),p4s=symbol("p4s"),p12s=symbol("p12s"),p23s=symbol("p23s");
	
	
	int c;
	int digit_optind = 0;
	double temp_numcast = 0;
	short int loopToolInput_Indicator = 0; // Whether input in P (as for loop tool):
	short int xloopsInput_Indicator = 0; // Whether input in Q (xloop type)
	
	
	while (1)
		{
			int this_option_optind = optind ? optind : 1;
			int option_index = 0;
			static struct option long_options[] =
			{
				{"m1s", 1, 0, 0},
				{"m2s", 1, 0, 0},
				{"m3s", 1, 0, 0},
				{"m4s", 1, 0, 0},
				{"p1s", 1, 0, 0},//P1 Square
				{"p2s", 1, 0, 0},//P2 Square
				{"p3s", 1, 0, 0},//P3 Square
				{"p4s", 1, 0, 0},//P4 Square
				{"p12s", 1, 0, 0},//(P1+P2) Square
				{"p23s", 1, 0, 0},//(P2+P3) Square
				{0, 0, 0, 0}
			};
			c = getopt_long (argc, argv, "",
				         long_options, &option_index);
			if (c == -1)
				break;
		
			switch (c)
			{
			case 0:
				temp_numcast = strtod(optarg, NULL);
				switch(option_index)
				{
				case 0: m1 = temp_numcast;
					break;
				case 1: m2 = temp_numcast;
					break;
				case 2: m3 = temp_numcast;
					break;
				case 3: m4 = temp_numcast;
					break;
				case 4: p1s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				case 5: p2s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				case 6: p3s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				case 7: p4s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				case 8: p12s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				case 9: p23s = temp_numcast; 
					loopToolInput_Indicator = 1;
					break;
				}
				break;
			default:
				throw std::runtime_error("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nUnknow input option, please use:\n --m1s --m2s --m3s --m4s --q10 --q20 --q21 --q30 --q31 --q32\nOR\n--m1s --m2s --m3s --m4s --p1s --p2s --p3s --p4s --p12s --p23s\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
			}
		}
			//input m2s, m3s, m4s, m1s for xloop
		m = lst(m2, m3, m4, m1);
		p = lst(p1s, p2s, p3s, p4s, p12s, p23s);
	}
/* ************************************************************ THE END *****************************************************************/



